Vignette to reproduce the graph in ‘How Social Learning Shapes the Efficacy of Preventative Health Behaviors in an Outbreak’

Main Paper

Fig2

A sigmoid function sig:

\[ p \sim \frac{1}{1+e^{-\kappa\times(x-\nu)}} \]

par(mar=c(4,4,1,1))
x=seq(0,1,.01)
inp=seq(0,1,.05)
clrs=colorRampPalette(c("purple4","yellow"))(length(inp))
plot(x,sig(x),type="n",ylim=c(0,1),xlim=c(0,1),xlab="% infected",ylab=expression(P[switch]),main="" )
for(i in 1:length(inp)) lines(x,sig(x,b=inp[i],a=4),ylim=c(0,1),xlim=c(0,1),col=clrs[i],lwd=2)
leg=seq(1,length(inp),length.out=5)
legend("bottomright",legend=rev(sapply(inp[leg],function(d)as.expression(bquote(.(d))))),title=expression(nu),col=rev(clrs[leg]),lwd=2,cex=.8,bg="white")

tstp=25
stp=rev(seq(-3,3,length.out=tstp))
clrs=colorRampPalette(c("purple4","yellow"))(tstp)
plot(x,sig(x),type="n",ylim=c(0,1),xlim=c(0,1),xlab="% infected",ylab=expression(P[switch]),main="" )
for(i in rev(1:length(stp))) lines(x,sig(x,a=10^stp[i]),ylim=c(0,1),xlim=c(0,1),col=clrs[i],lwd=2)
leg=seq(1,tstp,length.out=5)
legend("bottomright",legend=sapply(round(stp[leg]),function(d)as.expression(bquote(.(10^d)))),col=clrs[leg],lwd=2,title=expression(kappa),cex=.8,bg="white")

illustration of switch parameters

Fig3.tiff

A figure to illustrate different way to flatten the curve. We ran several simulationi in two different setup to show the differences between two situation: one when there is no social distancing, one in which there is.

Static simulations

The simulations have been run using:

library(parallel)
xsize=ysize=100
pop=generatePopulation(N=500,xsize=xsize,ysize=ysize,recovery=c(8,14)*25,speed=c(1,.2),behavior=rep(B,500))
cl <- makeForkCluster(18,outfile="")
neutralNSD=parSapply(cl,1:500,function(i){
                     print(i);
                     simu=slsirSimu(tstep=1500,p=c(1,1),di=2,i0=1,xsize=xsize,ysize=ysize,
                                    pop=pop,
                                    inf=9,
                                    sat=1000,
                                    p_i=1,
                                    strategy="best",
                                    ts=T,ap=F,visu=F
                                    )$timeseries[,2]

})
save(file="neutralNSD.bin",neutralNSD)

pop[,"behavior"]=G
neutralSD=parSapply(cl,1:500,function(i){
                    print(i);
                    simu=slsirSimu(tstep=1500,p=c(.1,.1),di=2,i0=1,xsize=xsize,ysize=ysize,
                                   pop=pop,
                                   inf=9,
                                   sat=1000,
                                   p_i=1,
                                   strategy="best",
                                   ts=T,ap=F,visu=F
                                   )$timeseries[,2]

})
save(file="neutralSD.bin",neutralSD)


neutralSD2=parSapply(cl,1:500,function(i){
                     print(i);
                     simu=slsirSimu(tstep=1500,p=c(.1,.1),di=2,i0=1,xsize=xsize,ysize=ysize,
                                    pop=pop,
                                    inf=9,
                                    sat=1000,
                                    p_i=1,
                                    strategy="random",
                                    ts=T,ap=F,visu=F
                                    )$timeseries[,2]

})
save(file="neutralSD2.bin",neutralSD2)
stopCluster(cl)

plot figure

par(mar=c(2,2,2,2))
plot(1,1,type="n",ylim=c(0,500),xlim=c(1,1500),ylab="",xlab="",xaxt="n",yaxt="n")
legend("topright",legend=c("No social distancing","Social distancing"),fill=c(myred,mygreen),cex=.95)
mtext("time",1,1)
mtext("number of infected people",2,1)

load("../data/neutralSD.bin")
load("../data/neutralNSD.bin")
load("../data/neutralSD2.bin")
#we use precomputed data to compute a mean curve, those data have been generated using the scripts in "exec/graphPaper.R"

meanNSD=apply(neutralNSD,1,mean)
meanSD=apply(neutralSD,1,mean)

polygon(0:1500,meanSD,col=adjustcolor(mygreen,.5),border=NA)
polygon(0:1500,meanNSD,col=adjustcolor(myred,.5),border=NA)
lines(meanSD,col=adjustcolor(mygreen,.9),lwd=2)
lines(meanNSD,col=adjustcolor(myred,.9),lwd=2)

abline(h=max(meanNSD),lty=2,col=myred)
abline(v=which.max(meanNSD),lty=2,col=myred)
abline(h=max(meanSD),lty=2,col=mygreen)
abline(v=which.max(meanSD),lty=2,col=mygreen)

text(bquote(paste("max infected (",I[max],") no SD")),x=1500,y=max(meanNSD)+10,pos=2,cex=.8)
mtext(bquote(paste("time max (",tau,") no SD")),at=which.max(meanNSD),adj=1,cex=.8)
text(bquote(paste("max infected (",I[max],") SD")),x=1500,y=max(meanSD)+10,pos=2,cex=.8)
mtext(bquote(paste("time max (",tau,") SD")),at=which.max(meanSD),adj=0,cex=.8)

abline(v=which.max(meanNSD),lty=2,col=myred)
abline(h=max(meanSD),lty=2,col=mygreen)
abline(v=which.max(meanSD),lty=2,col=mygreen)

Fig4.tiff

In this figure we present the overall results of the main simulations, with some some analysis of \(\delta\), the score of each simulation.

Run all simulations

To runn all simulation one need the script in exec

Those generate a lot of data, we then amde available the result of our simulation here:

you cann download

wget  httpe://adsdas.ads

To get the results of the simulation then just:

name='../data/midCurveAllBadBestSLS/'
aa=lapply(unlist(lapply(list.dirs(name),list.files,pattern="allresults.bin",full.names=T)),function(f){load(f);return(allresults)})
allresults=do.call("rbind",aa)
#allresults=allresults[allresults$inf <.02 & allresults$sat >10,]
allresults=allresults[-which(allresults$max_infect <4),] #we remove the results where no outbreak happend
allresults$distances=(1-bdistance(allresults$time_max) + bdistance(allresults$max_infect))/2

For this Figure we also want to compare how simulations with different score relates to the case with no social learning. For that we re-run simulations. These simulations take time so we give the results in data/

Re-run best/worst results

library(parallel)
cl <- makeForkCluster(16,outfile="")

best=allresults[order(allresults$distances),][1:1000,]
worst=allresults[order(allresults$distances,decreasing=T),][1:1000,]

repetBest=parSapply(cl,1:500,function(i){
                    j=sample(nrow(best),1);print(paste(i,j));
                    simu=slsirSimu(500,1500,p=c(1,.1),di=2,i0=1,recovery=c(8,14)*25,speed=c(1,.2),xsize=xsize,ysize=ysize,
                                inf=best$inf[j],
                                sat=best$sat[j],
                                inf_r=best$inf_r[j],
                                sat_r=best$sat_r[j],
                                p_i=best$pind[j],
                                sl_rad=best$sl_rad[j],
                                strategy="best",
                                ts=T,ap=F,visu=F
                                )$timeseries[,2]
})
save(file="repetBest.bin",repetBest)

repetWorst=parSapply(cl,1:500,function(i){
                    j=sample(nrow(worst),1);print(paste(i,j));
                    simu=slsirSimu(500,1500,p=c(1,.1),di=2,i0=1,recovery=c(8,14)*25,speed=c(1,.2),xsize=xsize,ysize=ysize,
                                inf=worst$inf[j],
                                sat=worst$sat[j],
                                inf_r=worst$inf_r[j],
                                sat_r=worst$sat_r[j],
                                p_i=worst$pind[j],
                                sl_rad=worst$sl_rad[j],
                                strategy="best",
                                ts=T,ap=F,visu=F
                                )$timeseries[,2]
})
save(file="repetWorst.bin",repetWorst)

stopCluster(cl)

load("../data/repetBest.bin") 
load("../data/repetWorst.bin") 
meanWorst=apply(repetWorst,1,mean) #
meanBest=apply(repetBest,1,mean) #

generate Figure

We can now plot the graph of the paper

### left panel
par(mar=c(4,4,2,2))
plot(1,1,type="n",ylim=c(0,500),xlim=c(1,1500),ylab="",xlab="",xaxt="n",yaxt="n")
abline(v=150,lwd=2,col="gray")
par(xpd=NA)
text(150,-20,"t=150",col="1",pos=1)
par(xpd=F)
legend("topright",legend=c("No Learning (all NA)","No Learning (all A)","Worst Simulations","Best Simulations"),col=c(myred,mygreen),lty=c(1,1,2,2),lwd=2,cex=.95)
axis(1)
axis(2)
mtext("time",1,3)
mtext("number of infected people",2,3)


lines(meanSD,col=adjustcolor(mygreen,.9),lwd=2)
lines(meanNSD,col=adjustcolor(myred,.9),lwd=2)
lines(meanBest,col=adjustcolor(mygreen,.9),lwd=2,lty=2)
lines(meanWorst,col=adjustcolor(myred,.9),lwd=2,lty=2)

### right panel
subresults=allresults[sample(nrow(allresults),20000),]
rank=rank(subresults$distances)
orderscol=rank
color_class=rev(brewer.pal(5,"PuBu"))
colorbest=color_class[1]
orderscol=rep(adjustcolor(color_class[5],.4),length(orderscol))
nvl=c(1000,2500,5000,10000)
for(i in (length(nvl):1)) orderscol[rank<=nvl[i]]=adjustcolor(color_class[i],.4)
#plot(subresults$time_max,subresults$max_infect,bg=orderscol,main=expression(I[max]*" wrt "*tau),pch=21,xlab=expression(tau),ylab=expression(I[max]),lwd=.1)
plot(subresults$time_max,subresults$max_infect,bg=orderscol,main="",pch=21,xlab=expression(tau),ylab=expression(I[max]),lwd=.1)
legend("topright",legend=c(paste("<",nvl),"all"),pt.bg=color_class,pch=21,col=adjustcolor(1,.5),pt.lwd=.1,title="Rank")

Fig 4: mean trajectories and scores for simulations

Fig5

For this figure we need the posterior distribution of our model given by the distance to the data as measured by our \(\delta\)

Get the posterior

n=1000
posterior=list(
               best=allresults[order(allresults$distances),][1:n,],
               worst=allresults[order(allresults$distances,decreasing=T),][1:n,],
               time_max=allresults[order(allresults$time_max,decreasing=T),][1:n,],
               time_max150=allresults[order(allresults$time_max150,decreasing=T),][1:n,],
               time_max250=allresults[order(allresults$time_max250,decreasing=T),][1:n,],
               wtime_max150=allresults[order(allresults$time_max150,decreasing=F),][1:n,],
               wtime_max250=allresults[order(allresults$time_max250,decreasing=F),][1:n,],
               max_infect=allresults[order(allresults$max_infect,decreasing=F),][1:n,],
               max_infect150=allresults[order(allresults$max_infect150,decreasing=F),][1:n,],
               max_infect250=allresults[order(allresults$max_infect250,decreasing=F),][1:n,],
               wmax_infect150=allresults[order(allresults$max_infect150,decreasing=T),][1:n,],
               wmax_infect250=allresults[order(allresults$max_infect250,decreasing=T),][1:n,]
               )

library(hdrcde)

duocolrs=adjustcolor(c(myred,colorbest),.8)
currange=range=list(dim1=c(0,1),dim2=c(-1,3))
dimlab=list(dim1=expression(nu),dim2=expression(kappa))

### Cmpound marginal with 2d  learning
expnames=c(distribA="Worst 150",distribB="Best 150")
marginAndJoin(
              distribA=list(dim1=posterior$wmax_infect150$inf,dim2=posterior$wmax_infect150$sat),
              distribB=list(dim1=posterior$max_infect150$inf,dim2=posterior$max_infect150$sat),
              dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
              expnames=expnames, range=currange,cols=duocolrs,log="y",
              main=""
              )

expnames=c(distribA="Worst 150",distribB="Best 150")
marginAndJoin(
              distribA=list(dim1=posterior$wmax_infect150$inf_r,dim2=posterior$wmax_infect150$sat_r),
              distribB=list(dim1=posterior$max_infect150$inf_r,dim2=posterior$max_infect150$sat_r),
              dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
              expnames=expnames, range=currange,cols=duocolrs,log="y",
              main=""
              )

expnames=c(distribA="Worst",distribB="Best")
marginAndJoin(
              distribA=list(dim1=posterior$worst$inf,dim2=posterior$worst$sat),
              distribB=list(dim1=posterior$best$inf,dim2=posterior$best$sat),
              dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
              expnames=expnames, range=currange,cols=duocolrs,log="y",
              main=""
              )


expnames=c(distribA="Worst",distribB="Best")
marginAndJoin(
              distribA=list(dim1=posterior$worst$inf_r,dim2=posterior$worst$sat_r),
              distribB=list(dim1=posterior$best$inf_r,dim2=posterior$best$sat_r),
              dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
              expnames=expnames, range=currange,cols=duocolrs,log="y",
              main=""
              )

Top row: distriution of parameter that gives us the best results when simulation reach 150 steps, bottom row: when the simulation reaches 1500. Left column: parameters to swtich from NA to A, right column to switch from NA to A

Supplementary Material

The code for the supplementary material generate multiples graphs, some have not been added into the paper.

Load data from all model:

##### Grpah papers
##get data
library(RColorBrewer)
modelnames=c("midCurveAllBad","midCurveAllBadBestSLS","midCurve15Good","midCurve15GoodBestSLS","midCurve100GoodBestSLS","burnin100BestSLS")
names(modelnames)=modelnames
modelcol=brewer.pal(length(modelnames),name="Dark2")
names(modelcol)=modelnames
modelsubnames=sapply(modelnames,function(mod)paste0(strsplit(mod,"")[[1]][-(1:8)],collapse=""))
modelsubnames=LETTERS[1:6]


#name="testradius"
setwd("../data/")

allexpe=lapply(modelnames,function(name){
               aa=lapply(unlist(lapply(list.dirs(name),list.files,pattern="allresults.bin",full.names=T)),function(f){load(f);return(allresults)})
               allresults=do.call("rbind",aa)
               if(sum(which(allresults$max_infect <4)>2))allresults=allresults[-which(allresults$max_infect <4),]
               allresults$distances=(1-bdistance(allresults$time_max) + bdistance(allresults$max_infect))/2
               allresults$distances2=(1-bdistance(allresults$time_max150) + bdistance(allresults$max_infect150))/2
               allresults$distances3=(1-bdistance(allresults$time_max250) + bdistance(allresults$max_infect250))/2
               return(allresults)})
nsims=sapply(allexpe,nrow)

#Bayes factors
allscores=c("distances","distances2","distances3","time_max","time_max150","time_max250","max_infect","max_infect150","max_infect250")
names(allscores)=allscores
tablesscores=sapply(allscores,function(as){
                    allscores=lapply(allexpe,"[[",as)
                    allnames=rep(names(allscores),lengths(allscores))
                    table(allnames[order(unlist(allscores),decreasing=F)][1:(10000)])
})
normalized=apply(tablesscores,2,"/",nsims[rownames(tablesscores)])
dis=normalized[,1]
##calculate bayes factor, by row, uie mat_rat[i,j] = dis[i]/dis[j]
mat_rat=sapply(1:length(dis),function(i)sapply(1:length(dis),function(j)dis[i]/dis[j]))  
colnames(mat_rat)=names(dis)
rownames(mat_rat)=names(dis)

Show distances for all models

lapply(names(allexpe),function(mname){
       allresults=allexpe[[mname]]
       plot(allresults$time_max,allresults$max_infect,col=adjustcolor(color.gradient(allresults$distances,c(mygreen,"yellow", myred)),.2),main=mname,pch=20,xlab="Time Max",ylab="Max Infected")
       sset=sort(allresults$distances)
       cols=adjustcolor(color.gradient(sset,c(mygreen,"yellow", myred)),1)
       lsset=seq(1,length(sset),length.out=4)
       legend("topright",legend=sapply(round(sset[lsset],digit=2),function(d)as.expression(bquote(delta==.(d)))),col=cols[lsset],pch=20)
})

## [[1]]
## [[1]]$rect
## [[1]]$rect$w
## [1] 140.4326
## 
## [[1]]$rect$h
## [1] 116.2025
## 
## [[1]]$rect$left
## [1] 883.3274
## 
## [[1]]$rect$top
## [1] 507.6
## 
## 
## [[1]]$text
## [[1]]$text$x
## [1] 926.5836 926.5836 926.5836 926.5836
## 
## [[1]]$text$y
## [1] 484.3595 461.1190 437.8785 414.6380
## 
## 
## 
## [[2]]
## [[2]]$rect
## [[2]]$rect$w
## [1] 151.7549
## 
## [[2]]$rect$h
## [1] 116.2025
## 
## [[2]]$rect$left
## [1] 936.4851
## 
## [[2]]$rect$top
## [1] 508.6
## 
## 
## [[2]]$text
## [[2]]$text$x
## [1] 983.2289 983.2289 983.2289 983.2289
## 
## [[2]]$text$y
## [1] 485.3595 462.1190 438.8785 415.6380
## 
## 
## 
## [[3]]
## [[3]]$rect
## [[3]]$rect$w
## [1] 143.5371
## 
## [[3]]$rect$h
## [1] 115.1772
## 
## [[3]]$rect$left
## [1] 900.9029
## 
## [[3]]$rect$top
## [1] 507.48
## 
## 
## [[3]]$text
## [[3]]$text$x
## [1] 945.1154 945.1154 945.1154 945.1154
## 
## [[3]]$text$y
## [1] 484.4446 461.4091 438.3737 415.3382
## 
## 
## 
## [[4]]
## [[4]]$rect
## [[4]]$rect$w
## [1] 146.0938
## 
## [[4]]$rect$h
## [1] 113.8101
## 
## [[4]]$rect$left
## [1] 909.9062
## 
## [[4]]$rect$top
## [1] 507.32
## 
## 
## [[4]]$text
## [[4]]$text$x
## [1] 954.9062 954.9062 954.9062 954.9062
## 
## [[4]]$text$y
## [1] 484.5580 461.7959 439.0339 416.2719
## 
## 
## 
## [[5]]
## [[5]]$rect
## [[5]]$rect$w
## [1] 140.7979
## 
## [[5]]$rect$h
## [1] 115.1772
## 
## [[5]]$rect$left
## [1] 886.0421
## 
## [[5]]$rect$top
## [1] 508.48
## 
## 
## [[5]]$text
## [[5]]$text$x
## [1] 929.4109 929.4109 929.4109 929.4109
## 
## [[5]]$text$y
## [1] 485.4446 462.4091 439.3737 416.3382
## 
## 
## 
## [[6]]
## [[6]]$rect
## [[6]]$rect$w
## [1] 154.8594
## 
## [[6]]$rect$h
## [1] 110.7342
## 
## [[6]]$rect$left
## [1] 978.0606
## 
## [[6]]$rect$top
## [1] 501.96
## 
## 
## [[6]]$text
## [[6]]$text$x
## [1] 1025.761 1025.761 1025.761 1025.761
## 
## [[6]]$text$y
## [1] 479.8132 457.6663 435.5195 413.3727

## distribution distances/metrcis/4models
par(oma=c(0,2,0,0),mar=c(4,1,1,0))
for(v in c(1,4,7)){
    var=allscores[v]
        xlimits=range(sapply(allexpe,"[[",var))
        alld=lapply(allexpe,function(i,var)density(i[[var]]),var=var)
        ylimits=range(sapply(alld,"[[","y"))
        t=""
        if(v==1)t=expression(delta)
        if(v==4)t=expression(tau)
        if(v==7)t=expression(I[max])
        plot(1,1,ylim=ylimits,xlim=xlimits,type="n",main=t,yaxt="n",xlab=t,ylab="density")
        lapply(names(alld),function(i)lines(alld[[i]],col=modelcol[i],lwd=2))

}
legend("topleft",legend=modelsubnames,col=modelcol,lwd=3,lty=1)
mtext("density",2,0,outer=T)

for(var in allscores[c(7)+1]){
        xlimits=range(sapply(allexpe,"[[",var))
        alld=lapply(allexpe,function(i,var)density(i[[var]]),var=var)
        ylimits=range(sapply(alld,"[[","y"))
        plot(1,1,ylim=ylimits,xlim=xlimits,type="n",main=var,yaxt="n",xlab=var,ylab="density")
        lapply(names(alld),function(i)lines(alld[[i]],col=modelcol[i],lwd=3))

}
legend("topright",legend=modelsubnames,col=modelcol,lwd=3,lty=1)

for(var in allscores[c(1,4,7)+2]){
        xlimits=range(sapply(allexpe,"[[",var))
        alld=lapply(allexpe,function(i,var)density(i[[var]]),var=var)
        ylimits=range(sapply(alld,"[[","y"))
        plot(1,1,ylim=ylimits,xlim=xlimits,type="n",main=var,yaxt="n",xlab=var,ylab="density")
        lapply(names(alld),function(i)lines(alld[[i]],col=modelcol[i],lwd=3))

}
legend("topright",legend=modelsubnames,col=modelcol,lwd=3,lty=1)

sapply(allexpe,function(i)quantile(i[[allscores[4]]],prob=.5)) 

##         midCurveAllBad.50%  midCurveAllBadBestSLS.50% 
##                      387.0                      384.5 
##         midCurve15Good.50%  midCurve15GoodBestSLS.50% 
##                      389.0                      385.0 
## midCurve100GoodBestSLS.50%       burnin100BestSLS.50% 
##                      387.0                      424.0

Posterior distribution for all models

n=1000
listallposterior=lapply(allexpe,function(e)
                        {
                            list(
                                 best=e[order(e$distances),][1:n,],
                                 best=e[order(e$distances),][1:n,],
                                 worst=e[order(e$distances,decreasing=T),][1:n,],
                                 time_max=e[order(e$time_max,decreasing=T),][1:n,],
                                 time_max150=e[order(e$time_max150,decreasing=T),][1:n,],
                                 time_max250=e[order(e$time_max250,decreasing=T),][1:n,],
                                 max_infect=e[order(e$max_infect,decreasing=F),][1:n,],
                                 max_infect150=e[order(e$max_infect150,decreasing=F),][1:n,],
                                 max_infect250=e[order(e$max_infect250,decreasing=F),][1:n,],
                                 wtime_max=e[order(e$time_max,decreasing=F),][1:n,],
                                 wtime_max150=e[order(e$time_max150,decreasing=F),][1:n,],
                                 wtime_max250=e[order(e$time_max250,decreasing=F),][1:n,],
                                 wmax_infect=e[order(e$max_infect,decreasing=T),][1:n,],
                                 wmax_infect150=e[order(e$max_infect150,decreasing=T),][1:n,],
                                 wmax_infect250=e[order(e$max_infect250,decreasing=T),][1:n,]
                                 )
                        }
)

Show the shape of sigmoid for the posteriors distributions

for(mname in names(listallposterior)){
    x=seq(0,1,.01)
    par(mfrow=c(length(listallposterior[[mname]]),2),oma=c(3,3,2,0))
    par(mar=rep(.5,4))
    for(mod in names(listallposterior[[mname]])){
        best=listallposterior[[mname]][[mod]]
        plot(x,sig(x),type="n",ylim=c(0,1),xlim=c(0,1),xlab="% same age infected people", ylab="proba. to switch",main="",xaxt="n",yaxt="n")
        modsubname=paste0(strsplit(mname,"")[[1]][-(1:8)],collapse="")
        axis(2)
        mtext(modsubname,2.5,2,cex=.9)
        a=axis(1,outer=T)
        for(i in 1:nrow(best)) lines(x,sig(x,b=best$inf[i],a=best$sat[i]),ylim=c(0,1),xlim=c(0,1),col=adjustcolor(modelcol[mname],.1),lwd=2) 
        plot(x,sig(x),type="n",ylim=c(0,1),xlim=c(0,1),xlab="% same age infected people", ylab="proba. to switch",main="",xaxt="n",yaxt="n")
        for(i in 1:nrow(best)) lines(x,sig(x,b=best$inf_r[i],a=best$sat_r[i]),ylim=c(0,1),xlim=c(0,1),col=adjustcolor(modelcol[mname],.1),lwd=2) 
    }
    axis(1,outer=F)
    mtext(mname,3,1,outer=T)
}

FIGURE 7

modsubname=sapply(modelnames,function(mod)paste0(strsplit(mod,"")[[1]][-(1:8)],collapse=""))
parameters=colnames(allexpe$midCurve15Good)[1:6]
par(mfrow=c(3,length(parameters)),mar=c(1,2,2,1))
for(best in listallposterior[6:8]){
    for(var in parameters){
        p=unlist(sapply(allexpe,"[[",var))
        xlimits=range(sapply(best,"[[",var),p)
        alld=lapply(best,function(i,var)density(i[[var]],from=xlimits[1],to=xlimits[2]),var=var)
        pd=density(p)
        if(var %in% c("sat_r","sat")){
            print(var)
            xlimits=c(-1,3)
            alld=lapply(best,function(i,var)density(log10(i[[var]]),from=xlimits[1],to=xlimits[2]),var=var)
            pd=density(log10(p),from=xlimits[1],to=xlimits[2])
        }
        ylimits=range(c(sapply(alld,"[[","y"),pd$y))
        plot(1,1,ylim=ylimits,xlim=xlimits,type="n",main=var,yaxt="n",xlab=var,ylab="density")
        lines(pd,lwd=3)
        lapply(names(alld),function(i)lines(alld[[i]],col=modelcol[i],lwd=3))
        if(var=="inf_r")legend("topleft",legend=modsubname,col=modelcol,lwd=3,lty=1)

    }
}

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

#### better separation use those::
parameters=colnames(allexpe$midCurve15Good)[1:6]
for(n in names(listallposterior)){
    sbest=listallposterior[[n]]
    #png(paste0("allposterior_metric_",n,".png"),width=1200,height=480,pointsize=15)
    par(mfrow=c(2,length(parameters)),mar=c(1,2,2,1),oma=c(3,3,0,0))
    sub=list(1:2,3:6)
    for( s in sub){
        best=sbest[s]
        for(var in parameters){
            p=unlist(sapply(allexpe,"[[",var))
            xlimits=range(sapply(best,"[[",var),p)
            alld=lapply(best,function(i,var)density(i[[var]],from=xlimits[1],to=xlimits[2]),var=var)
            pd=density(p)
            if(var %in% c("sat_r","sat")){
                print(var)
                xlimits=c(-1,3)
                alld=lapply(best,function(i,var)density(log10(i[[var]]),from=xlimits[1],to=xlimits[2]),var=var)
                pd=density(log10(p),from=xlimits[1],to=xlimits[2])
            }
            ylimits=range(c(sapply(alld,"[[","y"),pd$y))
            plot(1,1,ylim=ylimits,xlim=xlimits,type="n",main=var,yaxt="n",xlab=var,ylab="density")
            lines(pd,lwd=3)
            lapply(names(alld),function(i)lines(alld[[i]],col=modelcol[i],lwd=3))
            if(var=="sl_rad")legend("bottom",legend=modelsubnames[s],col=modelcol[s],lwd=3,lty=1)
        }

        }
        mtext("t0 : 15-100% A",2,1,0.3,outer=T) 
        mtext("t0 : 0% NA",2,1,0.7,outer=T) 
}

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

## [1] "sat"

## [1] "sat_r"

FIGURE 8

for(n in names(listallposterior)){
    sbest=listallposterior[[n]]
    sbest=sbest[-5]

    png(paste0("posterior2d_",n,"_inf_r_sat_r.png"),pointsize=15, res = 100, width = 6, height = 6, units = "in")
    par(mar=c(1,1,1,1),oma=c(3,3,3,3),xpd=T)
    par(mfrow=c(3,2))
    for(mod in names(sbest)){
        best=sbest[[mod]]
        if(mod=="burnin100BestSLS")plot.new()
        hdr.boxplot.2d(best$inf_r,log10(best$sat_r),prob=seq(20,100,10),shadecols=adjustcolor(modelcol[mod],1),xlim=c(0,1),ylim=c(-1,3),xlab="revert inf. point",ylab="revert steepness (log10)")
    }
    mtext("inf_r",1,2,outer=T)
    mtext("sat_r",2,2,outer=T)

    mtext("random SL",3,1,.3,outer=T)
    mtext("best SL",3,1,.7,outer=T)

    mtext("wait100A",4,1,.15,outer=T)
    mtext("15%A",4,1,.5,outer=T)
    mtext("0%AD",4,1,,.85,outer=T)
    dev.off()

    png(paste0("posterior2d_",n,"_inf_sat.png"),pointsize=15, res = 100, width = 6, height = 6, units = "in")
    par(mar=c(1,1,1,1),oma=c(3,3,3,3),xpd=T)
    par(mfrow=c(3,2))
    for(mod in names(sbest)){
        best=listallposterior$allbest[[mod]]
        if(mod=="burnin100BestSLS")plot.new()
        hdr.boxplot.2d(best$inf,log10(best$sat),prob=seq(20,100,10),shadecols=adjustcolor(modelcol[mod],.9),xlim=c(0,1),ylim=c(-1,3),xlab="inflexion point",ylab="steepness (log10)")
    }
    mtext("inf",1,2,outer=T)
    mtext("sat",2,2,outer=T)
                                     
    mtext("random SL",3,1,.3,outer=T)
    mtext("best SL",3,1,.7,outer=T)

    mtext("wait100A",4,1,.15,outer=T)
    mtext("15%A",4,1,.5,outer=T)
    mtext("0%AD",4,1,,.85,outer=T)
                                 
    dev.off()

    png(paste0("posterior2d_",n,"_inf_pind.png"),pointsize=15, res = 100, width = 6, height = 6, units = "in")
    par(mar=c(1,1,1,1),oma=c(3,3,3,3),xpd=T)
    par(mfrow=c(3,2))
    for(mod in names(sbest)){
        best=listallposterior$allbest[[mod]]
        if(mod=="burnin100BestSLS")plot.new()
        hdr.boxplot.2d(best$inf,best$pind,prob=seq(20,100,10),shadecols=adjustcolor(modelcol[mod],.9),xlim=c(0,1),ylim=c(0,1),xlab="inflexion point",ylab="individual learning")
    }
    mtext("inf",1,2,outer=T)
    mtext("pind",2,2,outer=T)

    mtext("random SL",3,1,.3,outer=T)
    mtext("best SL",3,1,.7,outer=T)

    mtext("wait100A",4,1,.15,outer=T)
    mtext("15%A",4,1,.5,outer=T)
    mtext("0%AD",4,1,,.85,outer=T)
    dev.off()

    png(paste0("posterior2d_",n,"_sat_pind.png"),pointsize=15, res = 100, width = 6, height = 6, units = "in")

    par(mar=c(1,1,1,1),oma=c(3,3,3,3),xpd=T)
    par(mfrow=c(3,2))
    for(mod in names(sbest)){
        best=listallposterior$allbest[[mod]]
        if(mod=="burnin100BestSLS")plot.new()
        hdr.boxplot.2d(log10(best$sat),best$pind,prob=seq(20,100,10),shadecols=adjustcolor(modelcol[mod],.9),xlim=c(-1,3),ylim=c(0,1),xlab="steepness (log10)",ylab="individual learning")
    }
    mtext("sat",1,2,outer=T)
    mtext("pind",2,2,outer=T)

    mtext("random SL",3,1,.3,outer=T)
    mtext("best SL",3,1,.7,outer=T)
                                
    mtext("wait100A",4,1,.15,outer=T)
    mtext("15%A",4,1,.5,outer=T)
    mtext("0%AD",4,1,,.85,outer=T)
                                 
    dev.off()

    png(paste0("posterior2d_",n,"_sl_rad_pind.png"),pointsize=15, res = 100, width = 6, height = 6, units = "in")
    par(mar=c(1,1,1,1),oma=c(3,3,3,3),xpd=T)
    par(mfrow=c(3,2))
    for(mod in names(sbest)){
        best=listallposterior$allbest[[mod]]
        if(mod=="burnin100BestSLS")plot.new()
        hdr.boxplot.2d(best$sl_rad,best$pind,prob=seq(20,100,10),shadecols=adjustcolor(modelcol[mod],.9),xlim=c(1,100),ylim=c(0,1),xlab="radius social learning",ylab="individual learning")
    }
    mtext("sl_rad",1,2,outer=T)
    mtext("pind",2,2,outer=T)
    mtext("random SL",3,1,.3,outer=T)
    mtext("best SL",3,1,.7,outer=T)
    
    mtext("wait100A",4,1,.15,outer=T)
    mtext("15%A",4,1,.5,outer=T)
    mtext("0%AD",4,1,,.85,outer=T)
    dev.off()

    png(paste0("posterior2d_",n,"_sl_rad_inf.png"),pointsize=15, res = 100, width = 6, height = 6, units = "in")
    par(mar=c(1,1,1,1),oma=c(3,3,3,3),xpd=T)
    par(mfrow=c(3,2))
    for(mod in names(sbest)){
        best=listallposterior$allbest[[mod]]
        if(mod=="burnin100BestSLS")plot.new()
        hdr.boxplot.2d(best$sl_rad,best$inf,prob=seq(20,100,10),shadecols=adjustcolor(modelcol[mod],.9),xlim=c(1,100),ylim=c(0,1),xlab="radius social learning",ylab="infection points")
    }
    mtext("sl_rad",1,2,outer=T)
    mtext("inf",2,2,outer=T)

    mtext("random SL",3,1,.3,outer=T)
    mtext("best SL",3,1,.7,outer=T)
                                
    mtext("wait100A",4,1,.15,outer=T)
    mtext("15%A",4,1,.5,outer=T)
    mtext("0%AD",4,1,,.85,outer=T)
                                 
    dev.off()

    png(paste0("posterior2d_",n,"_sl_rad_sat.png"),pointsize=15, res = 100, width = 6, height = 6, units = "in")
par(xpd=NA)
    par(mar=c(1,1,1,1),oma=c(3,3,3,3),xpd=T)
    par(mfrow=c(3,2))
    for(mod in names(sbest)){
        best=listallposterior$allbest[[mod]]
        if(mod=="burnin100BestSLS")plot.new()
        hdr.boxplot.2d(best$sl_rad,log10(best$sat),prob=seq(20,100,10),shadecols=adjustcolor(modelcol[mod],.9),xlim=c(1,100),ylim=c(-1,3),xlab="sl_rad",ylab="steepness (log10)")
    }
    mtext("sl_rad",1,2,outer=T)
    mtext("sat",2,2,outer=T)

    mtext("random SL",3,1,.3,outer=T)
    mtext("best SL",3,1,.7,outer=T)
                                
    mtext("wait100A",4,1,.15,outer=T)
    mtext("15%A",4,1,.5,outer=T)
    mtext("0%AD",4,1,,.85,outer=T)
                                 
    dev.off()
}

par(mfrow=c(2,4),oma=c(0,4,0,0))
 
plot(density(log10(listallposterior$allbest$midCurveAllBad$sat),from=-1,to=3),ylim=c(0,.45),main="sat",xlab="sat")
lines(density(log10(listallposterior$allbest_max_infect150$midCurveAllBad$sat)),col="red")
lines(density(log10(allexpe$midCurveAllBad$sat)),col="grey")
legend("topright",legend=c("final time","150ts", "prior"),lwd=1,col=c(1,2,"grey"))
mtext("normal model",2,5)

plot(density(listallposterior$allbest$midCurveAllBad$inf,from=0,to=1),ylim=c(0,2),main="inf",xlab="inf")
lines(density(listallposterior$allbest_max_infect150$midCurveAllBad$inf,from=0,to=1),col="red")
lines(density(allexpe$midCurveAllBad$inf),col="grey")

plot(density(log10(listallposterior$allbest$midCurveAllBad$sat_r),from=-1,to=3),ylim=c(0,.7),main="sat_r",xlab="sat_r")
lines(density(log10(listallposterior$allbest_max_infect150$midCurveAllBad$sat_r)),col="red")
lines(density(log10(allexpe$midCurveAllBad$sat_r)),col="grey")

plot(density(listallposterior$allbest$midCurveAllBad$inf_r,from=0,to=1),ylim=c(0,6),main="inf_r",xlab="inf_r")
lines(density(listallposterior$allbest_max_infect150$midCurveAllBad$inf_r,from=0,to=1),col="red")
lines(density(allexpe$midCurveAllBad$inf_r),col="grey")

plot(density(log10(listallposterior$allbest$midCurveAllBad$sat),from=-1,to=3),ylim=c(0,.45),main="sat",xlab="sat")
lines(density(log10(listallposterior$allbest_max_infect150$midCurveAllBad$sat)),col="red")
lines(density(log10(allexpe$midCurveAllBad$sat)),col="grey")
legend("topright",legend=c("final time","150ts", "prior"),lwd=1,col=c(1,2,"grey"))
mtext("waitn 100ts ",2,5)

plot(density(listallposterior$allbest$midCurveAllBad$inf,from=0,to=1),ylim=c(0,2),main="inf",xlab="inf")
lines(density(listallposterior$allbest_max_infect150$midCurveAllBad$inf,from=0,to=1),col="red")
lines(density(allexpe$midCurveAllBad$inf),col="grey")

plot(density(log10(listallposterior$allbest$midCurveAllBad$sat_r),from=-1,to=3),ylim=c(0,.7),main="sat_r",xlab="sat_r")
lines(density(log10(listallposterior$allbest_max_infect150$midCurveAllBad$sat_r)),col="red")
lines(density(log10(allexpe$midCurveAllBad$sat_r)),col="grey")

plot(density(listallposterior$allbest$midCurveAllBad$inf_r,from=0,to=1),ylim=c(0,6),main="inf_r",xlab="inf_r")
lines(density(listallposterior$allbest_max_infect150$midCurveAllBad$inf_r,from=0,to=1),col="red")
lines(density(allexpe$midCurveAllBad$inf_r),col="grey")



par(mfrow=c(1,2))
plot(1,1,ylim=c(0,500),xlim=c(0,500),type="n",main="#infected")
 for(f in list.files("testingSimHigSL/" ,pattern="allG_*",full.names=T)){ 
     load(f)
     lines(simu$timeseries[,2],col="green")
}
 for(f in list.files("testingSimHigSL/" ,pattern="allB_*",full.names=T)){ 
     load(f)
     lines(simu$timeseries[,2],col="red")
}
plot(1,1,ylim=c(0,500),xlim=c(0,500),type="n",main="#Adherent")
 for(f in list.files("testingSimHigSL/" ,pattern="allG_*",full.names=T)){ 
     load(f)
     lines(simu$timeseries[,5],col="green")
}
 for(f in list.files("testingSimHigSL/" ,pattern="allB_*",full.names=T)){ 
     load(f)
     lines(simu$timeseries[,5],col="red")
}

Compound marginal with 2d learning

duocolrs=adjustcolor(c(myred,colorbest),.8)
currange=range=list(dim1=c(0,1),dim2=c(-1,3))
dimlab=list(dim1=expression(nu),dim2=expression(kappa))

for(i in names(allexpe)){
    posterior=listallposterior[[i]]
    worst=allexpe[[i]]

    expnames=c(distribA="Worst",distribB="Best")
    marginAndJoin(
                  distribA=list(dim1=posterior$worst$inf,dim2=posterior$worst$sat),
                  distribB=list(dim1=posterior$best$inf,dim2=posterior$best$sat),
                  dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )

    expnames=c(distribA="Worst 150",distribB="Best 150")
    marginAndJoin(
                  distribA=list(dim1=posterior$wmax_infect150$inf,dim2=posterior$wmax_infect150$sat),
                  distribB=list(dim1=posterior$max_infect150$inf,dim2=posterior$max_infect150$sat),
                  dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )

    ### COmpound marginal with 2d revert learning

    expnames=c(distribA="Worst",distribB="Best")
    marginAndJoin(
                  distribA=list(dim1=posterior$worst$inf_r,dim2=posterior$worst$sat_r),
                  distribB=list(dim1=posterior$best$inf_r,dim2=posterior$best$sat_r),
                  dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )

    expnames=c(distribA="Worst 150",distribB="Best 150")
    marginAndJoin(
                  distribA=list(dim1=posterior$wmax_infect150$inf_r,dim2=posterior$wmax_infect150$sat_r),
                  distribB=list(dim1=posterior$max_infect150$inf_r,dim2=posterior$max_infect150$sat_r),
                  dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )

    ################ best vs prior

    ### Cmpound marginal with 2d  learning

    expnames=c(distribA="Prior",distribB="Best")
    marginAndJoin(
                  distribA=list(dim1=worst$inf,dim2=worst$sat),
                  distribB=list(dim1=posterior$best$inf,dim2=posterior$best$sat),
                  dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )

    expnames=c(distribA="Prior",distribB="Best 150")
    marginAndJoin(
                  distribA=list(dim1=worst$inf,dim2=worst$sat),
                  distribB=list(dim1=posterior$max_infect150$inf,dim2=posterior$max_infect150$sat),
                  dimlab=list(dim1=expression(nu),dim2=expression(kappa)),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )

    ### COmpound marginal with 2d revert learning

    expnames=c(distribA="Prior",distribB="Best")
    marginAndJoin(
                  distribA=list(dim1=worst$inf_r,dim2=worst$sat_r),
                  distribB=list(dim1=posterior$best$inf_r,dim2=posterior$best$sat_r),
                  dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )

    expnames=c(distribA="Prior",distribB="Best 150")
    marginAndJoin(
                  distribA=list(dim1=worst$inf_r,dim2=worst$sat_r),
                  distribB=list(dim1=posterior$max_infect150$inf_r,dim2=posterior$max_infect150$sat_r),
                  dimlab=list(dim1=expression(nu[r]),dim2=expression(kappa[r])),
                  expnames=expnames, range=currange,cols=duocolrs,log="y",
                  main=""
                  )
}

FIGURE N

for(i in names(allexpe)){
    posterior=listallposterior[[i]]
    mar=c(4,3,2,1)

    colors=c(P="NA",A=adjustcolor(colorbest,1),B=adjustcolor(colorbest,1))

    par(mar=mar)
    plot2dens(prior=1-allresults$pind,A=1-posterior$best$pind,B=1-posterior$max_infect250$pind,cols=colors,from=0,to=1,main="posterior proba social learning",xlab="proba social learning")

    plot2dens(prior=allresults$sl_rad,A=posterior$best$sl_rad,B=posterior$max_infect150$sl_rad,cols=colors,from=1,to=100,main="",xlab="radius social learning")

    plot2dens(prior=log10(allresults$sat_r),A=log10(posterior$best$sat_r),B=log10(posterior$time_max150$sat_r),cols=colors,from=-1,to=3,main="",xlab=expression(kappa[r]),log=T)

    plot2dens(prior=log10(allresults$sat),A=log10(posterior$best$sat),B=log10(posterior$max_infect150$sat),cols=colors,from=-1,to=3,xlab=expression(kappa),main="",xaxt="n",log=T)

    plot2dens(prior=(allresults$inf_r),A=(posterior$best$inf_r),B=(posterior$max_infect150$inf_r),cols=colors,from=0,to=1,main="posterior revert inflection point",xlab=expression(nu[r]),xaxt="n")
    legend("topleft",legend=c("prior","final time","150 timesteps"),fill=c(0,colorbest,colorbest),density=c(NA,NA,20))
    
    plot2dens(prior=(allresults$inf),A=(posterior$best$inf),B=(posterior$max_infect150$inf),cols=colors,from=0,to=1,main="posterior inflection point",xlab=expression(nu),xaxt="n")
}